// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2008 Gael Guennebaud <gael.guennebaud@inria.fr>
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.

#ifndef EIGEN_AMBIVECTOR_H
#define EIGEN_AMBIVECTOR_H

namespace Eigen {

namespace internal {

    /** \internal
  * Hybrid sparse/dense vector class designed for intensive read-write operations.
  *
  * See BasicSparseLLT and SparseProduct for usage examples.
  */
    template <typename _Scalar, typename _StorageIndex> class AmbiVector
    {
    public:
        typedef _Scalar Scalar;
        typedef _StorageIndex StorageIndex;
        typedef typename NumTraits<Scalar>::Real RealScalar;

        explicit AmbiVector(Index size) : m_buffer(0), m_zero(0), m_size(0), m_end(0), m_allocatedSize(0), m_allocatedElements(0), m_mode(-1) { resize(size); }

        void init(double estimatedDensity);
        void init(int mode);

        Index nonZeros() const;

        /** Specifies a sub-vector to work on */
        void setBounds(Index start, Index end)
        {
            m_start = convert_index(start);
            m_end = convert_index(end);
        }

        void setZero();

        void restart();
        Scalar& coeffRef(Index i);
        Scalar& coeff(Index i);

        class Iterator;

        ~AmbiVector() { delete[] m_buffer; }

        void resize(Index size)
        {
            if (m_allocatedSize < size)
                reallocate(size);
            m_size = convert_index(size);
        }

        StorageIndex size() const { return m_size; }

    protected:
        StorageIndex convert_index(Index idx) { return internal::convert_index<StorageIndex>(idx); }

        void reallocate(Index size)
        {
            // if the size of the matrix is not too large, let's allocate a bit more than needed such
            // that we can handle dense vector even in sparse mode.
            delete[] m_buffer;
            if (size < 1000)
            {
                Index allocSize = (size * sizeof(ListEl) + sizeof(Scalar) - 1) / sizeof(Scalar);
                m_allocatedElements = convert_index((allocSize * sizeof(Scalar)) / sizeof(ListEl));
                m_buffer = new Scalar[allocSize];
            }
            else
            {
                m_allocatedElements = convert_index((size * sizeof(Scalar)) / sizeof(ListEl));
                m_buffer = new Scalar[size];
            }
            m_size = convert_index(size);
            m_start = 0;
            m_end = m_size;
        }

        void reallocateSparse()
        {
            Index copyElements = m_allocatedElements;
            m_allocatedElements = (std::min)(StorageIndex(m_allocatedElements * 1.5), m_size);
            Index allocSize = m_allocatedElements * sizeof(ListEl);
            allocSize = (allocSize + sizeof(Scalar) - 1) / sizeof(Scalar);
            Scalar* newBuffer = new Scalar[allocSize];
            std::memcpy(newBuffer, m_buffer, copyElements * sizeof(ListEl));
            delete[] m_buffer;
            m_buffer = newBuffer;
        }

    protected:
        // element type of the linked list
        struct ListEl
        {
            StorageIndex next;
            StorageIndex index;
            Scalar value;
        };

        // used to store data in both mode
        Scalar* m_buffer;
        Scalar m_zero;
        StorageIndex m_size;
        StorageIndex m_start;
        StorageIndex m_end;
        StorageIndex m_allocatedSize;
        StorageIndex m_allocatedElements;
        StorageIndex m_mode;

        // linked list mode
        StorageIndex m_llStart;
        StorageIndex m_llCurrent;
        StorageIndex m_llSize;
    };

    /** \returns the number of non zeros in the current sub vector */
    template <typename _Scalar, typename _StorageIndex> Index AmbiVector<_Scalar, _StorageIndex>::nonZeros() const
    {
        if (m_mode == IsSparse)
            return m_llSize;
        else
            return m_end - m_start;
    }

    template <typename _Scalar, typename _StorageIndex> void AmbiVector<_Scalar, _StorageIndex>::init(double estimatedDensity)
    {
        if (estimatedDensity > 0.1)
            init(IsDense);
        else
            init(IsSparse);
    }

    template <typename _Scalar, typename _StorageIndex> void AmbiVector<_Scalar, _StorageIndex>::init(int mode)
    {
        m_mode = mode;
        // This is only necessary in sparse mode, but we set these unconditionally to avoid some maybe-uninitialized warnings
        // if (m_mode==IsSparse)
        {
            m_llSize = 0;
            m_llStart = -1;
        }
    }

    /** Must be called whenever we might perform a write access
  * with an index smaller than the previous one.
  *
  * Don't worry, this function is extremely cheap.
  */
    template <typename _Scalar, typename _StorageIndex> void AmbiVector<_Scalar, _StorageIndex>::restart() { m_llCurrent = m_llStart; }

    /** Set all coefficients of current subvector to zero */
    template <typename _Scalar, typename _StorageIndex> void AmbiVector<_Scalar, _StorageIndex>::setZero()
    {
        if (m_mode == IsDense)
        {
            for (Index i = m_start; i < m_end; ++i) m_buffer[i] = Scalar(0);
        }
        else
        {
            eigen_assert(m_mode == IsSparse);
            m_llSize = 0;
            m_llStart = -1;
        }
    }

    template <typename _Scalar, typename _StorageIndex> _Scalar& AmbiVector<_Scalar, _StorageIndex>::coeffRef(Index i)
    {
        if (m_mode == IsDense)
            return m_buffer[i];
        else
        {
            ListEl* EIGEN_RESTRICT llElements = reinterpret_cast<ListEl*>(m_buffer);
            // TODO factorize the following code to reduce code generation
            eigen_assert(m_mode == IsSparse);
            if (m_llSize == 0)
            {
                // this is the first element
                m_llStart = 0;
                m_llCurrent = 0;
                ++m_llSize;
                llElements[0].value = Scalar(0);
                llElements[0].index = convert_index(i);
                llElements[0].next = -1;
                return llElements[0].value;
            }
            else if (i < llElements[m_llStart].index)
            {
                // this is going to be the new first element of the list
                ListEl& el = llElements[m_llSize];
                el.value = Scalar(0);
                el.index = convert_index(i);
                el.next = m_llStart;
                m_llStart = m_llSize;
                ++m_llSize;
                m_llCurrent = m_llStart;
                return el.value;
            }
            else
            {
                StorageIndex nextel = llElements[m_llCurrent].next;
                eigen_assert(i >= llElements[m_llCurrent].index && "you must call restart() before inserting an element with lower or equal index");
                while (nextel >= 0 && llElements[nextel].index <= i)
                {
                    m_llCurrent = nextel;
                    nextel = llElements[nextel].next;
                }

                if (llElements[m_llCurrent].index == i)
                {
                    // the coefficient already exists and we found it !
                    return llElements[m_llCurrent].value;
                }
                else
                {
                    if (m_llSize >= m_allocatedElements)
                    {
                        reallocateSparse();
                        llElements = reinterpret_cast<ListEl*>(m_buffer);
                    }
                    eigen_internal_assert(m_llSize < m_allocatedElements && "internal error: overflow in sparse mode");
                    // let's insert a new coefficient
                    ListEl& el = llElements[m_llSize];
                    el.value = Scalar(0);
                    el.index = convert_index(i);
                    el.next = llElements[m_llCurrent].next;
                    llElements[m_llCurrent].next = m_llSize;
                    ++m_llSize;
                    return el.value;
                }
            }
        }
    }

    template <typename _Scalar, typename _StorageIndex> _Scalar& AmbiVector<_Scalar, _StorageIndex>::coeff(Index i)
    {
        if (m_mode == IsDense)
            return m_buffer[i];
        else
        {
            ListEl* EIGEN_RESTRICT llElements = reinterpret_cast<ListEl*>(m_buffer);
            eigen_assert(m_mode == IsSparse);
            if ((m_llSize == 0) || (i < llElements[m_llStart].index))
            {
                return m_zero;
            }
            else
            {
                Index elid = m_llStart;
                while (elid >= 0 && llElements[elid].index < i) elid = llElements[elid].next;

                if (llElements[elid].index == i)
                    return llElements[m_llCurrent].value;
                else
                    return m_zero;
            }
        }
    }

    /** Iterator over the nonzero coefficients */
    template <typename _Scalar, typename _StorageIndex> class AmbiVector<_Scalar, _StorageIndex>::Iterator
    {
    public:
        typedef _Scalar Scalar;
        typedef typename NumTraits<Scalar>::Real RealScalar;

        /** Default constructor
      * \param vec the vector on which we iterate
      * \param epsilon the minimal value used to prune zero coefficients.
      * In practice, all coefficients having a magnitude smaller than \a epsilon
      * are skipped.
      */
        explicit Iterator(const AmbiVector& vec, const RealScalar& epsilon = 0) : m_vector(vec)
        {
            using std::abs;
            m_epsilon = epsilon;
            m_isDense = m_vector.m_mode == IsDense;
            if (m_isDense)
            {
                m_currentEl = 0;    // this is to avoid a compilation warning
                m_cachedValue = 0;  // this is to avoid a compilation warning
                m_cachedIndex = m_vector.m_start - 1;
                ++(*this);
            }
            else
            {
                ListEl* EIGEN_RESTRICT llElements = reinterpret_cast<ListEl*>(m_vector.m_buffer);
                m_currentEl = m_vector.m_llStart;
                while (m_currentEl >= 0 && abs(llElements[m_currentEl].value) <= m_epsilon) m_currentEl = llElements[m_currentEl].next;
                if (m_currentEl < 0)
                {
                    m_cachedValue = 0;  // this is to avoid a compilation warning
                    m_cachedIndex = -1;
                }
                else
                {
                    m_cachedIndex = llElements[m_currentEl].index;
                    m_cachedValue = llElements[m_currentEl].value;
                }
            }
        }

        StorageIndex index() const { return m_cachedIndex; }
        Scalar value() const { return m_cachedValue; }

        operator bool() const { return m_cachedIndex >= 0; }

        Iterator& operator++()
        {
            using std::abs;
            if (m_isDense)
            {
                do
                {
                    ++m_cachedIndex;
                } while (m_cachedIndex < m_vector.m_end && abs(m_vector.m_buffer[m_cachedIndex]) <= m_epsilon);
                if (m_cachedIndex < m_vector.m_end)
                    m_cachedValue = m_vector.m_buffer[m_cachedIndex];
                else
                    m_cachedIndex = -1;
            }
            else
            {
                ListEl* EIGEN_RESTRICT llElements = reinterpret_cast<ListEl*>(m_vector.m_buffer);
                do
                {
                    m_currentEl = llElements[m_currentEl].next;
                } while (m_currentEl >= 0 && abs(llElements[m_currentEl].value) <= m_epsilon);
                if (m_currentEl < 0)
                {
                    m_cachedIndex = -1;
                }
                else
                {
                    m_cachedIndex = llElements[m_currentEl].index;
                    m_cachedValue = llElements[m_currentEl].value;
                }
            }
            return *this;
        }

    protected:
        const AmbiVector& m_vector;  // the target vector
        StorageIndex m_currentEl;    // the current element in sparse/linked-list mode
        RealScalar m_epsilon;        // epsilon used to prune zero coefficients
        StorageIndex m_cachedIndex;  // current coordinate
        Scalar m_cachedValue;        // current value
        bool m_isDense;              // mode of the vector
    };

}  // end namespace internal

}  // end namespace Eigen

#endif  // EIGEN_AMBIVECTOR_H
